library(TSA)
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(fpp)
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.Arima TSA
## plot.Arima TSA
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: tseries
library(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(readxl)
library(corrplot)
## corrplot 0.92 loaded
library(stats)
df_raw <- read_excel('TS_Project.xlsx')
df <- df_raw
head(df)
## # A tibble: 6 × 9
## time export import cpi ppi bond sentiment USDX uncertai…¹
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2002-01-01 00:00:00 1.41 3.57 178. 128. 5.04 93 120. 117.
## 2 2002-02-01 00:00:00 1.34 3.52 178 128. 4.91 90.7 119. 87.6
## 3 2002-03-01 00:00:00 1.76 3.79 178. 130. 5.28 95.7 119. 83.4
## 4 2002-04-01 00:00:00 1.46 4.45 179. 131. 5.21 93 115. 88.6
## 5 2002-05-01 00:00:00 1.63 4.36 180. 131. 5.16 96.9 112. 86.3
## 6 2002-06-01 00:00:00 1.76 4.48 180. 131. 4.93 92.4 106. 93.6
## # … with abbreviated variable name ¹​uncertainty
tail(df)
## # A tibble: 6 × 9
## time export import cpi ppi bond sentiment USDX uncertai…¹
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2022-09-01 00:00:00 6.72 28.5 297. 268. 3.52 58.6 112. 174.
## 2 2022-10-01 00:00:00 5.91 28.4 298. 265. 3.98 59.9 112. 177.
## 3 2022-11-01 00:00:00 5.50 25.6 299. 263. 3.89 56.8 106. 172.
## 4 2022-12-01 00:00:00 5.85 24.8 299. 258. 3.62 59.7 104. 136.
## 5 2023-01-01 00:00:00 6.16 25.0 301. 260. 3.53 64.9 102. 143.
## 6 2023-02-01 00:00:00 5.86 20.6 302. 259. 3.75 67 105. 124.
## # … with abbreviated variable name ¹​uncertainty
any(is.na(df))
## [1] FALSE
export = ts(df$export, start = c(2002, 1), frequency = 12)
import = ts(df$import, start = c(2002, 1), frequency = 12)
cpi = ts(df$cpi, start=c(2002,1),frequency=12)
ppi = ts(df$ppi, start=c(2002,1),frequency=12)
bond = ts(df$bond, start=c(2002,1),frequency=12)
sentiment = ts(df$sentiment, start=c(2002,1),frequency=12)
usdx = ts(df$USDX, start=c(2002,1),frequency=12)
uncertainty = ts(df$uncertainty, start=c(2002,1),frequency=12)
# create a correlation matrix
cor_matrix <- cor(cbind(export, import,
cpi, ppi,
bond, sentiment,
usdx, uncertainty))
# Plot a correlation heatmap
corrplot(cor_matrix, method = "color",
order = "hclust",
tl.col = "black", tl.srt = 45, addCoef.col = "black")
# Print the correlation table
# CPI & PPI highly correlated, only keep one.
print(cor_matrix)
## export import cpi ppi bond
## export 1.0000000 0.9576931 0.9508205 0.88258986 -0.61290296
## import 0.9576931 1.0000000 0.9623065 0.91101757 -0.64355565
## cpi 0.9508205 0.9623065 1.0000000 0.92833707 -0.69498684
## ppi 0.8825899 0.9110176 0.9283371 1.00000000 -0.56578495
## bond -0.6129030 -0.6435556 -0.6949868 -0.56578495 1.00000000
## sentiment -0.1779160 -0.2435011 -0.2186952 -0.43137804 0.07714775
## usdx 0.2673713 0.2693136 0.2465156 0.01599541 -0.07202055
## uncertainty 0.3651635 0.4453334 0.4491599 0.38277243 -0.59587323
## sentiment usdx uncertainty
## export -0.17791601 0.26737131 0.36516347
## import -0.24350112 0.26931356 0.44533338
## cpi -0.21869519 0.24651557 0.44915994
## ppi -0.43137804 0.01599541 0.38277243
## bond 0.07714775 -0.07202055 -0.59587323
## sentiment 1.00000000 0.32526534 -0.46620569
## usdx 0.32526534 1.00000000 0.04882897
## uncertainty -0.46620569 0.04882897 1.00000000
autoplot(import, main="Import") + ylab('Import(in billion $)')
autoplot(cpi, main="CPI") + ylab('CPI')
autoplot(ppi, main="PPI") + ylab('PPI')
autoplot(bond, main="10Y Govn't Bond") + ylab('Bond')
autoplot(sentiment, main='Consumer Sentiment') + ylab('Sentiment')
autoplot(usdx, main='USD Index') + ylab('USD Index')
autoplot(uncertainty, main='Economic Uncertainty Index') + ylab('Economic Uncertainty Index')
lambda <- BoxCox.lambda(import)
lambda
## [1] -0.08103269
cbind('Before Transformation' = import,
'After Transformation'= BoxCox(import, lambda = lambda)) %>% autoplot(facet = TRUE) + ggtitle('BoxCox Transformation Comparison') + ylab('Import(B$)')
# KPSS Test with the original data
kpss.test(import)
## Warning in kpss.test(import): p-value smaller than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: import
## KPSS Level = 3.7518, Truncation lag parameter = 5, p-value = 0.01
# KPSS Test with 1st order differencing data
diff1 = diff(import)
kpss.test(diff1)
## Warning in kpss.test(diff1): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff1
## KPSS Level = 0.027579, Truncation lag parameter = 5, p-value = 0.1
plot(diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing")
# KPSS Test with 1st order differncing + baxcox data
bc_diff1 = diff(BoxCox(import, lambda = lambda))
kpss.test(bc_diff1)
## Warning in kpss.test(bc_diff1): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: bc_diff1
## KPSS Level = 0.10007, Truncation lag parameter = 5, p-value = 0.1
plot(bc_diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing with BoxCox Transformation")
# Original Data
Acf(import)
Pacf(import)
# Data with 1st order differencing
Acf(diff1)
Pacf(diff1)
# Data with 1st order differencig + boxcox
Acf(bc_diff1)
Pacf(bc_diff1)
# Original Plot
autoplot(import)
# Decomposition Plot (STL)
ts_data <- ts(bc_diff1, frequency = 12)
stl_fit <- stl(ts_data, s.window = "periodic")
plot(stl_fit)
remainder <- stl_fit$time.series[, "remainder"]
plot(remainder)
# Seasonal Plot
ggseasonplot(bc_diff1, polar = TRUE)
ggseasonplot(bc_diff1, polar = FALSE)
# Subseries Plot
ggsubseriesplot(bc_diff1) +
labs(y = "Import(B$)",
title = "Seasonal plot: Import")
# Periodogram
spectrum(bc_diff1)
periodogram(bc_diff1)
temp <- periodogram(bc_diff1)
max_freq <- temp$freq[which.max(temp$spec)]
seasonality <- 1/max_freq
1 / temp$freq[43] # 5.953488
## [1] 5.953488
1 / temp$freq[107] #2.392523
## [1] 2.392523
1 / temp$freq[21] # 12.19048
## [1] 12.19048
Acf(bc_diff1)
Pacf(bc_diff1)
df_raw <- read_excel('TS_Project.xlsx')
df <- ts(df_raw$import, start=c(2002,01), frequency=12)
autoplot(df) + ylab('Import') + ggtitle('US Monthly Import(2002/01 ~ 2023/02)')
df_train <- window(df, start=c(2002, 1), end=c(2018, 12))
df_test <- window(df, start=c(2019, 1), end=c(2019, 12))
fcst_snaive <- snaive(df_train, h = 12)
autoplot(snaive(df_train, h = 12))
accuracy(fcst_snaive)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.7617909 1.527044 1.280903 7.1367 12.59435 1 0.7750122
accuracy(fcst_snaive$mean, df_test)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.4499564 1.248044 1.03557 2.415875 5.877225 -0.3960617 0.7297898
checkresiduals(fcst_snaive)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 682.64, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
ses3 <- ets(df_train, model='ANN')
fcst_ses3 <- forecast(ses3, h=12)
accuracy(fcst_ses3)[,2]
## [1] 0.9509481
accuracy(fcst_ses3$mean, df_test)[,2]
## [1] 1.643461
checkresiduals(fcst_ses3$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 172.5, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
autoplot(ses3) + ylab('Import($B)')
autoplot(ses3) + autolayer(fitted(ses3)) + ylab('Import($B)')
# 1) Linear trend with additive seasonality
hwa <- hw(df_train, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
## ACF1
## Training set 0.01406249
accuracy(hwa$mean, df_test)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.1271767 0.8469332 0.6890742 0.4701921 3.986005 -0.4069091 0.4849101
checkresiduals(hwa$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
##
## Model df: 0. Total lags used: 24
# pdf("plot.pdf", width = 10, height = 6)
plot(hwa)
lines(df_test, col = "red")
# 2) Linear trend with multiplicative seasonality
hwm <- hw(df_train, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02226477 0.5818108 0.4363622 -0.01005463 4.344781 0.3406676
## ACF1
## Training set -0.03132167
accuracy(hwm$mean, df_test)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.2187908 0.8660047 0.7331691 0.9624011 4.256856 -0.3067592 0.5090764
checkresiduals(hwm$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 42.431, df = 24, p-value = 0.01154
##
## Model df: 0. Total lags used: 24
plot(hwm)
lines(df_test, col = "red")
# 3) Linear trend with additive seasonality and damping
hwad <- hw(df_train, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.08943375 0.6597823 0.4971829 0.8109482 5.319108 0.3881503
## ACF1
## Training set -0.01316333
accuracy(hwad$mean, df_test)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.2356866 1.105793 0.913462 0.7802561 5.422637 0.1745505 0.6544695
checkresiduals(hwad$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.299, df = 24, p-value = 0.01986
##
## Model df: 0. Total lags used: 24
plot(hwad)
lines(df_test, col = "red")
# 4) Linear trend with multiplicative seasonality and damping
hwmd <- hw(df_train, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03776369 0.5807925 0.4377261 0.2958402 4.469252 0.3417324
## ACF1
## Training set -0.02054866
accuracy(hwmd$mean, df_test)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.4950831 0.9971667 0.8049686 2.566052 4.561863 -0.2181502 0.5892602
checkresiduals(hwmd$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 35.775, df = 24, p-value = 0.0577
##
## Model df: 0. Total lags used: 24
plot(hwmd)
lines(df_test, col = "red")
# 1 ETS_BEST
ets_best <- ets(df_train)
summary(ets_best)
## ETS(M,A,M)
##
## Call:
## ets(y = df_train)
##
## Smoothing parameters:
## alpha = 0.5353
## beta = 1e-04
## gamma = 2e-04
##
## Initial states:
## l = 4.0986
## b = 0.0424
## s = 0.9908 1.0682 1.1106 1.0167 1.0183 1.0399
## 1.0286 0.9756 0.9869 0.9801 0.8526 0.9316
##
## sigma: 0.057
##
## AIC AICc BIC
## 847.9923 851.2826 904.4003
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03831786 0.5906786 0.4378753 0.1381647 4.319084 0.3418489
## ACF1
## Training set -0.03521243
checkresiduals(ets_best$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 42.295, df = 24, p-value = 0.01196
##
## Model df: 0. Total lags used: 24
fcst_ets_best <- forecast(ets_best, h=12)
accuracy(fcst_ets_best$mean, df_test)[,2] #RMSE:2.534726
## [1] 0.8525379
autoplot(forecast(ets_best)) + ylab('Import($B)')
plot(ets_best)
plot(fcst_ets_best)
lines(df_test, col = "red")
# 2. ETS_AAA
ets1 <- ets(df_train, model='AAA')
summary(ets1)
## ETS(A,A,A)
##
## Call:
## ets(y = df_train, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.4797
## beta = 1e-04
## gamma = 0.2997
##
## Initial states:
## l = 4.0752
## b = 0.0591
## s = -0.2185 0.1891 0.3965 -0.0095 0.2708 0.3481
## 0.3534 0.1724 0.1472 -0.1551 -0.8346 -0.6597
##
## sigma: 0.6759
##
## AIC AICc BIC
## 942.3971 945.6874 998.8051
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
## ACF1
## Training set 0.01406249
checkresiduals(ets1$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
##
## Model df: 0. Total lags used: 24
fcst_ets1 <- forecast(ets1, h=12)
accuracy(fcst_ets1$mean, df_test)[,2]
## [1] 0.8469332
autoplot(forecast(ets1)) + ylab('Import($B)')
plot(ets1)
plot(fcst_ets1)
lines(df_test, col = "red")
# 3. ETS MAA
ets2 <- ets(df_train, model='MAA')
summary(ets2)
## ETS(M,A,A)
##
## Call:
## ets(y = df_train, model = "MAA")
##
## Smoothing parameters:
## alpha = 0.5215
## beta = 1e-04
## gamma = 0.206
##
## Initial states:
## l = 4.1148
## b = 0.0575
## s = 0.0136 0.3469 0.5217 0.1403 0.2663 0.2594
## 0.1876 -0.0763 -0.0638 -0.1533 -0.8209 -0.6217
##
## sigma: 0.0615
##
## AIC AICc BIC
## 880.9684 884.2588 937.3765
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.007855664 0.6566577 0.4792403 -0.2734595 4.698261 0.3741425
## ACF1
## Training set 0.004567237
checkresiduals(ets2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 38.875, df = 24, p-value = 0.02815
##
## Model df: 0. Total lags used: 24
fcst_ets2 <- forecast(ets2, h=12)
accuracy(fcst_ets2$mean, df_test)[,2]
## [1] 0.8176148
autoplot(forecast(ets2)) + ylab('Import($B)')
plot(ets2)
plot(fcst_ets2)
lines(df_test, col = "red")
# 4. ETS MMM
ets3 <- ets(df_train, model='MMM')
summary(ets3)
## ETS(M,M,M)
##
## Call:
## ets(y = df_train, model = "MMM")
##
## Smoothing parameters:
## alpha = 0.554
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 4.037
## b = 1.0072
## s = 0.9947 1.0738 1.1076 1.0177 1.0197 1.0366
## 1.0237 0.9757 0.9852 0.9811 0.8512 0.933
##
## sigma: 0.0567
##
## AIC AICc BIC
## 847.6910 850.9813 904.0990
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01799567 0.5900927 0.4405025 -0.2720982 4.360535 0.3438999
## ACF1
## Training set -0.05031675
checkresiduals(ets3)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,M,M)
## Q* = 40.548, df = 24, p-value = 0.01866
##
## Model df: 0. Total lags used: 24
fcst_ets3 <- forecast(ets3, h=12)
accuracy(fcst_ets3$mean, df_test)[,2]
## [1] 0.8680253
autoplot(forecast(ets3)) + ylab('Import($B)')
plot(ets3)
plot(fcst_ets3)
lines(df_test, col = "red")
df_train2 <- window(df, start=c(2002, 1), end=c(2021, 12))
df_test2 <- window(df, start=c(2022, 1), end=c(2022,12))
fcst_snaive <- snaive(df_train2, h = 12)
autoplot(snaive(df_train2, h = 12))
accuracy(fcst_snaive $mean, df_test2)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 3.242702 4.445163 4.111537 12.07819 15.53626 0.4679377 1.611265
ses32 <- ets(df_train2, model='ANN')
fcst_ses32 <- forecast(ses32, h=12)
accuracy(fcst_ses32)[,2]
## [1] 1.118179
accuracy(fcst_ses32$mean, df_test2)[,2]
## [1] 2.103524
checkresiduals(fcst_ses32$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 211.62, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
autoplot(ses32) + ylab('Import($B)')
autoplot(ses32) + autolayer(fitted(ses32)) + ylab('Import($B)')
# 1) Linear trend with additive seasonality
hwa2 <- hw(df_train2, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
## ACF1
## Training set 0.01213894
accuracy(hwa2$mean, df_test2)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.142104 3.432513 2.447793 -8.391731 9.549002 0.4576512 1.273803
checkresiduals(hwa2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
##
## Model df: 0. Total lags used: 24
plot(hwa2)
lines(df_test2, col = "red")
# 2) Linear trend with multiplicative seasonality
hwm2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04448036 0.7014787 0.5107636 0.005129007 4.569962 0.3460351
## ACF1
## Training set 0.06717551
accuracy(hwm2$mean, df_test2)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.570837 3.174066 2.495035 -6.166118 9.709162 0.4946886 1.172481
checkresiduals(hwm2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 60.478, df = 24, p-value = 5.467e-05
##
## Model df: 0. Total lags used: 24
plot(hwm2)
lines(df_test2, col = "red")
# 3) Linear trend with additive seasonality and damping
hwad2 <- hw(df_train2, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.07959986 0.7972801 0.6008245 0.5927588 5.847758 0.4070501
## ACF1
## Training set 0.05177696
accuracy(hwad2$mean, df_test2)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -3.092347 3.829269 3.108496 -12.07427 12.13198 0.2978517 1.418878
checkresiduals(hwad2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 58.148, df = 24, p-value = 0.0001161
##
## Model df: 0. Total lags used: 24
plot(hwad2)
lines(df_test2, col = "red")
# 4) Linear trend with multiplicative seasonality and damping
hwmd2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd2)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05315586 0.6975741 0.5102869 0.1244391 4.548272 0.3457122
## ACF1
## Training set 0.02301767
accuracy(hwmd2$mean, df_test2)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.176211 3.633945 2.760502 -8.450697 10.74315 0.5221199 1.339611
checkresiduals(hwmd2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 52.418, df = 24, p-value = 0.0006899
##
## Model df: 0. Total lags used: 24
plot(hwmd2)
lines(df_test2, col = "red")
# 1. ETS_BEST
ets_best2 <- ets(df_train2)
summary(ets_best2)
## ETS(M,A,M)
##
## Call:
## ets(y = df_train2)
##
## Smoothing parameters:
## alpha = 0.5311
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 4.0927
## b = 0.0597
## s = 1.0063 1.0846 1.1101 1.0194 1.0135 1.0378
## 1.0187 0.9762 0.9892 0.9744 0.8362 0.9336
##
## sigma: 0.0581
##
## AIC AICc BIC
## 1094.530 1097.287 1153.701
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06211406 0.7294272 0.5263934 0.0311818 4.498561 0.3566241
## ACF1
## Training set 0.01879104
checkresiduals(ets_best2$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 46.18, df = 24, p-value = 0.004212
##
## Model df: 0. Total lags used: 24
fcst_ets_best2 <- forecast(ets_best2, h=12)
accuracy(fcst_ets_best2$mean, df_test2)[,2]
## [1] 2.019115
autoplot(forecast(ets_best2)) + ylab('Import($B)')
plot(ets_best2)
plot(fcst_ets_best2)
lines(df_test2, col = "red")
# 2. ETS_AAA
ets12 <- ets(df_train2, model='AAA')
summary(ets12)
## ETS(A,A,A)
##
## Call:
## ets(y = df_train2, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.4867
## beta = 0.0276
## gamma = 0.346
##
## Initial states:
## l = 4.2205
## b = 0.1176
## s = 0.202 0.0634 0.5544 -0.1505 -0.2504 0.3038
## 0.2454 0.0457 0.3895 0.1847 -0.6942 -0.8938
##
## sigma: 0.7971
##
## AIC AICc BIC
## 1223.920 1226.677 1283.091
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
## ACF1
## Training set 0.01213894
checkresiduals(ets12$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
##
## Model df: 0. Total lags used: 24
fcst_ets12 <- forecast(ets12, h=12)
accuracy(fcst_ets12$mean, df_test2)[,2]
## [1] 3.432513
autoplot(forecast(ets12)) + ylab('Import($B)')
plot(ets12)
plot(fcst_ets12)
lines(df_test2, col = "red")
# 3. ETS MAA
ets22 <- ets(df_train2, model='MAA')
summary(ets22)
## ETS(M,A,A)
##
## Call:
## ets(y = df_train2, model = "MAA")
##
## Smoothing parameters:
## alpha = 0.5565
## beta = 1e-04
## gamma = 0.2409
##
## Initial states:
## l = 4.3895
## b = 0.0618
## s = -0.042 0.2136 0.6139 0.1723 0.2651 0.2256
## 0.2715 -0.1323 0.0125 -0.2621 -0.7543 -0.5837
##
## sigma: 0.062
##
## AIC AICc BIC
## 1126.354 1129.110 1185.524
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05355054 0.773105 0.5608316 -0.1563987 4.828425 0.3799555
## ACF1
## Training set 0.01192749
checkresiduals(ets22$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 35.481, df = 24, p-value = 0.06159
##
## Model df: 0. Total lags used: 24
fcst_ets22 <- forecast(ets22, h=12)
accuracy(fcst_ets22$mean, df_test2)[,2]
## [1] 1.907392
autoplot(forecast(ets22)) + ylab('Import($B)')
plot(ets22)
plot(fcst_ets22)
lines(df_test2, col = "red")
# 4. ETS MMM
ets32 <- ets(df_train2, model='MMM')
summary(ets32)
## ETS(M,M,M)
##
## Call:
## ets(y = df_train2, model = "MMM")
##
## Smoothing parameters:
## alpha = 0.5646
## beta = 1e-04
## gamma = 9e-04
##
## Initial states:
## l = 4.0127
## b = 1.0074
## s = 1.0144 1.0841 1.1101 1.0237 1.0135 1.0326
## 1.0172 0.9665 0.9824 0.9872 0.838 0.9304
##
## sigma: 0.058
##
## AIC AICc BIC
## 1094.544 1097.301 1153.715
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01542948 0.7196395 0.5242459 -0.1484539 4.482695 0.3551692
## ACF1
## Training set -0.03362053
checkresiduals(ets32$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 46.808, df = 24, p-value = 0.003535
##
## Model df: 0. Total lags used: 24
fcst_ets32 <- forecast(ets32, h=12)
accuracy(fcst_ets32$mean, df_test2)[,2]
## [1] 2.680248
autoplot(forecast(ets32)) + ylab('Import($B)')
plot(ets32)
plot(fcst_ets32)
lines(df_test2, col = "red")